Astron. Nachr. / AN 326, No. 9, (2005) / DOI 10. 1002/asna.2005 10414 



Turbulence and its parameterization in accretion discs 



Axel Brandenburg 

NORDITA, Blegdamsvej 17, DK-2100 Copenhagen 0, Denmaik 



Received 8 September 2005; accepted 22 September 2005; published online 20 October 2005 



Abstract. Accretion disc turbulence is investigated in the framework of the shearing box approximation. The turbulence is 
either driven by the magneto-rotational instability or, in the non-magnetic case, by an explicit and artificial forcing term in 
the momentum equation. Unlike the magnetic case, where most of the dissipation occurs in the disc corona, in the forced 
hydrodynamic case most of the dissipation occurs near the midplane. In the hydrodynamic case evidence is presented for the 
stochastic excitation of epicycles. When the vertical and radial epicyclic frequencies are different (modeling the properties 
around rotating black holes), the beat frequency between these two frequencies appear to show up as a peak in the temporal 
power spectrum in some cases. Finally, the full turbulent resistivity tensor is determined and it is found that, if the turbulence 
is driven by a forcing term, the signs of its off-diagonal components are such that this effect would not be capable of dynamo 
action by the shear-current effect. 
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1. Introduction 

It is now generally accepted that angular momentum transport 
in accretion discs is accomplished by hydromagnetic turbu- 
lence that is produced by the magneto-rotational instability 
(MRI, also known as Balbus-Hawley instability); see Balbus 
& Hawley (1991, 1998). The significance of this mechanism 
for accretion discs has been established using local shear- 
ing box simulations (Hawley et al. 1995, 1996, Matsumoto 
& Tajima 1995, Brandenburg et al. 1995, 1996a, Stone et 
al. 1996), as well as global simulations (Hawley 2000, Arlt 
& Rudiger 2001, De Villiers & Hawley 2003). For many 
purposes one would like to parameterize the turbulence in 
terms of a turbulent viscosity. The ultimate goal of such an 
approach is to be able to capture the relevant pieces of tur- 
bulence physics in two-dimensional axisymmetric and one- 
dimensional vertically integrated models of accretion discs. 
Even if this turns out not to be possible, parameterized mod- 
els are still extremely useful for illuminating otherwise unrec- 
ognized mechanisms that might only be directly identifiable 
using a targeted approach. 

One aspect that we do not understand very well right now 
is to what extent MRI-driven turbulence is similar to ordinary 
(e.g. forced) turbulence. This question is relevant because one 
is tempted to apply some well-known turbulence concepts 
quite loosely also to MRI-driven turbulence without distin- 



guishing between the different forms of driving. In order to 
address this question we consider here the case of forced tur- 
bulence and compare with what has been found for MRI- 
driven turbulence. In addition to the forced simulations we 
also present solutions of MRI-driven turbulence that have not 
previously been published. No net magnetic flux is imposed, 
so we are able to have a self-sustained mechanism. Most 
previous approaches used numerical viscosity and resistivity. 
With the regular laplacian diffusion operator (proportional to 
V^) self-sustained turbulence is only possible at large reso- 
lution (more than 128"^ meshpoints; see Brandenburg et al. 
2004). This becomes prohibitively expensive if we need to 
achieve sufficient statistics and long enough runs. Therefore, 
we adopt hyperdiffusion for these runs (here proportional to 
V^), analogous to what was done also in Brandenburg et al. 
(1995). All forced simulations are however done with regular 
laplacian diffusion. 

One of the goals of this paper is to reconsider the numer- 
ical determination of turbulent viscosity and resistivity in lo- 
cal simulations of accretion flows. This is motivated by recent 
advances in the case of non-shearing and non-rotating flows. 
Particular attention is paid to the tensorial nature of the turbu- 
lent resistivity tensor and the differences between forced and 
MRI-driven turbulence. For the turbulent viscosity the tenso- 
rial nature is important for modeling warps in accretion discs 
(Torkelsson et al. 2000). This has also motivated the study of 
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the tensorial nature of turbulent passive scalar diffusion (Car- 
ballido et al. 2005, Johansen & Klahr 2005). 

In the following we use overbars to denote spatial aver- 
ages over one or two coordinate directions, e.g. azimuthal 
averages and occasionally also vertical averages. Angular 
brackets without subscript are used for volume averages, 
while angular brackets with subscript t denote time averages. 

2. Stress and strain 

In any parameterized model of a turbulent flow one is inter- 
ested in the Reynolds stress 'pujuj and, if magnetic fields 

are present, in the Maxwell stress ^Sijb^ — bibj. Here, the 
magnetic field is measured in units where the vacuum per- 
meability is unity, and lower case characters u and b denote 
deviations from the mean flow U and the mean field B, so 
U = U + u and B = B + b m'e the full velocity and mag- 
netic fields. So the full turbulent stress from the small scales 
(SS) is given by 

^^(SS) 1 ^ J 2 1 — J — /I \ 

iiy = pUiUj + ^dijb - bibj. (1) 

Analogously, one can define the total stress from the large 
scale (LS) fields, 

n^^'^'^ = pUdJj + - B[Bj . (2) 

In the steady state, the value of the turbulent mass accretion 
rate, for example, follows from the constancy of the angular 
momentum flux, i.e. 

^ wd<j,j^dz (ntf + nL''^) = C = const, (3) 

where we have neglected the microscopic viscosity, because 
it is very small in discs (although not necessarily in simu- 
lations!). Here, cylindrical polar coordinates, (ro, 0, z), have 
been employed, and h denotes the disc height (e.g. its gaus- 
sian scale height). In Eq. (jSj we can isolate the mass accre- 
tion rate, M = — j vjdcf) j dzpUT^. Replacing furthermore 
the integration by a multiplication with A-Kwh, we find 

Af = (vo-ntf - voBjB^ - C) . (4) 

Remarkably, for calculating the accretion rate it suffices to 
know only the value of the stress, not its functional depen- 
dence on the mean flow properties. The same is true of the 
heating rate, which is given by (e.g., Balbus & Papaloizou 
1999, Balbus 2004) 

E^«.(.g)(nLT + nS>), (5, 

where 51 — /m has been introduced, and only the shear 
from the differential rotation has been taken into account. 

However, for many (if not all) other purposes it is neces- 
sary to know also the functional dependence of the stress on 
other quantities, most notably the shear rate. In fact, a com- 
mon proposal is to approximate the mean stress by the mean 
rate-of-strain tensor Sij and to write 

nir^ ~ 2pz^tSij, (6) 
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Fig. 1. Joule dissipation for Run 3 of Brandenburg (2001, 
solid line), compared with the Joule dissipation estimated for 
a corresponding mean-field model (dashed line). An estimate 
for the rate of total energy dissipation, B1^/t, is also given. 
[Adapted from Brandenburg (2003).] 

where Sy = \{Ui.j + Uj,i) is the rate of strain of the mean 
flow, and vt is a turbulent transport coefficient (turbulent vis- 
cosity).' 

The significance of Eq. (|6j is that it provides a closure 

(SS) 

for the small scale quantity 11^ in terms of the large scale 
strain, Sy . This is, unlike the previous relations in this sec- 
tion, necessarily only an approximation. In the equation of 
angular momentum conservation this term acts as a diffusion 
term and provides angular momentum transport in the direc- 
tion of decreasing angular velocity, i.e. outward. 

As we discussed above, the stress also contributes to the 
rate of heating of the disc. The horizontally averaged rate of 
viscous heating (per unit volume) is given by 

Qvisc = 2pj/S^ (actual heating rate), (7) 

where S is the actual rate of strain matrix, is assumed to be 
approximated by 

—2 

Qvisc ~ 2pj/tS (parameterized heating rate), (8) 

where S is the rate of strain of the mean flow. [In accretion 
discs, where keplerian shear gives the largest contribution to 
S , Eq. (|8} yields the familiar expression Qvisc ~ P^ti^^Y", 
see Frank et al. (1992).] Of course, » and ut ^ u, but 
whether vtS is actually the same as vS^ is not obvious. 

Using the first order smoothing approximation, which is 
commonly used in mean field dynamo theory, Riidiger (1987) 
found that vS^ is actually about 3 times larger than what is 

2 

expected from v^S .It is still unclear whether this is actually 
true, or whether it is an artifact of the first order smoothing 
approximation and that the two expressions should really be 
the same. 

In quite a different context of hydromagnetic turbulence, 
where a helical large scale mean field is generated (so that the 

' The mean flow is here solenoidal, so Sij is trace-less. In the 
general case considered below an extra term is added to make sure 
is trace-free. 
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Fig. 2. Dependence of the stress component (here denoted by r^y), separately for the kinetic and magnetic contri- 

butions, together with the sum of the two denoted by total (left) as well as the vertical dependence of density and sound 
speed (right). Note that Txy is neither proportional to the density p nor the sound speed Cs- [Adapted from Brandenburg et al. 
(1996b).] 



mean current density is well defined), an equivalent conclu- 
sion was reached for the actual and parameterized resistive 
heating rates, i.e. 



77J2 « 3.7 X ryt/, 



(9) 

see Brandenburg (2003) and Fig. where we also com- 
pare with an estimate for the rate of total energy dissipation, 
B^q/r, where r is the turnover time. Here we have used for 
rjt the value estimated from the self-consistently determined 
empirical quenching formula of Brandenburg (2001) for his 
Run 3. The factor 3.7 in Eq. (|9} is similar to what has been 
found for the turbulent viscosity (Riidiger 1987). This sug- 
gests that his result obtained within the framework of the first 
order smoothing approximation is at least not in conflict with 
the simulations. 

We now turn to another aspect of viscous heating. In 
MRI-driven turbulence it was found that the stress does not 
decrease with height away from the midplane, as suggested 
by Eq. (|6j, even though the product of rate of strain and 
density does decrease because of decreasing density. This 
was originally demonstrated only for nearly isothermal discs 
(Brandenburg et al. 1996b), see Fig.|2] but this has now also 
been shown for radiating discs (Turner 2004). 

Indeed, it was found that p should rather be replaced by 
the vertically averaged density. A sensible parameterization 
seemed therefore only possible for the vertically integrated 
stress, i.e. (Brandenburg et al. 1996b) 

d^^n^P « ass / dzpcg, (10) 

where ass is the dimensionless Shakura-Sunyaev (1973) vis- 
cosity coefficient and Cs is the sound speed. The reason for 
this is quite clearly related to the fact that the stress is of mag- 
netic origin, and that most of the dissipation happens when 
the density is low, i.e. when the ohmic heating rate per unit 
mass, j'^ / (crp), is large (here, J is the current density, and a 
the electric conductivity). However, it would be surprising if 
the absence of a scaling of the stress with the density were a 
general property. 
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Fig. 3. Vertical dependence of ass for hydrodynamic forced 
shear flow turbulence (upper panel). The solid curves give the 
time average while a few individual times are plotted in gray. 
Note that the temporal average of ass is nearly constant in 
z, although for individual times there can be significant fluc- 
tuations. By comparison, the density-weighted viscosity pa- 
rameter decreases with height (lower panel), showing that in 
forced non-magnetic turbulence the stress is indeed propor- 
tional to the density. 
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In order to shed some Hght on this question we now con- 
sider a hydrodynamic shearing box model where the tur- 
bulence is driven by an explicit body force that is delta- 
correlated in time and monochromatic in space with a 
wavenumber of k[ « 3. (Details of this simulation are given 
in Sect.|3]) The result is shown in Fig.|3l where we have ex- 
pressed the result in terms of a vertically dependent Shakura- 
Sunyaev viscosity coefficient ass(2)^ so 

vt{z) = assiz) Csh. (11) 

Note that in the present case of purely hydrodynamic forced 
turbulence ass{z) is nearly independent of height, even 
though p varies by at least an order of magnitude, as can be 
seen from the second panel of Fig. |3j where we show the 
product p{z)ass(z), normalized by the initial density in the 
midplane, po- 

In the following section we reiterate the standard shearing 
box equation that have been used to obtain this result. 

3. Shearing box equations 

We consider here a domain of size Lx x Ly x Lz, where 
Lx = Ly = 27T and Lz = 4, although other choices are 
possible and have been considered in the earlier work. We 
solve the equations for the departure from the keplerian shear 
flow for an isothermal equation of state (so the pressure is 
given by p — pc^). The resulting equations can be written in 
the form 



J X B 



Vlnp - 2fl X u 



Dt 

Dlnp 
Dt 



-A, 



dxi 



+ 



dxi 



3 —V■U+2S^J—-^ 

OXi OXi 



-V ■ tt, 



.(12) 



(13) 



(14) 



where D/Dt = d/dt+{U+u)-'V is the advective derivative 
based on the full flow field that includes both the shear flow 
and the deviations from it. We solve the equations for the 
departure, u, from the purely linear shear flow. Except for 
the advection operator, only derivatives of U enter This is an 
important property of the shearing sheet approximation that 
is critical for being able to use shearing-periodic boundary 
conditions in the x direction. We have used the gauge 

(j) — ryV • A — [U + m) ■ A, (assuming -q = const). (15) 

for the electrostatic potential. The shear flow is given by 
U{x) — (0, — qfiXjO), where q — 3/2 for a purely keple- 
rian shear flow, and 



duj 
dx-i 



dxi 



-5i,V ■ u 



(16) 



is the traceless rate of strain matrix. The Coriolis force is 
added to take into account that the shearing box is spinning 
about the central star at angular velocity fl. This together with 
shear can lead to radial epicyclic oscillations with frequency 

K= y/2(2^q)Q. (17) 



Note that k = fl for a keplerian disc with q = 3/2. The 
term — C^i; characterizes the vertical stratification, where ( is 
the vertical epicyclic frequency. Again, in a keplerian disc we 
have C = SI, but here we treat ( as an independent parame- 
ter in order to assess the effects of different radial and ver- 
tical epicyclic frequencies in non-newtonian discs with dif- 
ferent radial and vertical epicyclic frequencies (Abramowicz 
et al. 2003a,b, Kato 2004, Kluzniak et al. 2004, Lee et al. 
2004). For most of the cases considered below we choose 
n = K = 0.4 and ( = 0.6. The sound speed is taken to be 
Cs — 1- (With these parameters the magnitude of the shear 
flow is ^Q{Lx/2) Rs 1.9, so it is weakly supersonic.) In 
the nonmagnetic cases {A = 0) we allow for the possibil- 
ity of an extra forcing term / in order to study the case of 
non-magnetic (non-MRI) turbulence. The forcing function is 
identical to that used by Brandenburg (2001) for forced tur- 
bulent simulations exhibiting dynamo action. The amplitude 
of the forcing function is denoted by /o which is here chosen 
to be 0.01, unless it put to zero. The details of this forcing 
function are summarized in Appendix A. 

The boundary conditions adopted in the vertical direction 
on z = ±Lz/2 are 

"Ux.z = '^y.z = Uz = Ax = Ay = Az.z ~ 0, (18) 

which corresponds to a perfect conductor no-slip condition. 
Here, commas denote partial differentiation. 

The results presented here have been obtained with the 
Pencil Code,^ which is a public domain code for solv- 
ing the compressible hydromagnetic (and other) equations 
on massively parallel distributed memory architectures such 
as Beowulf clusters. It employs a high-order finite-difference 
scheme (sixth order in space and third order in time), which 
is ideal for all types of turbulence simulations. 

The present simulations are mainly exploratory in nature, 
so we only use small to moderate resolution between 32"^ and 
128^ meshpoints. Restrictions in resolution are also imposed 
by the need to run for long enough times in order to overcome 
transients. 

4. Growth and decay of epicycles 

One way of estimating the effective viscosity of a system is 
to determine the decay rate of a velocity perturbation that has 
initially a sinusoidal variation in one spatial direction. As an 
example one may consider an initial mean velocity perturba- 
tion of the form 

w(a;,0) = t7o(0)sinfc, -a; (19) 

for different directions i of the wavevector fe, . The time de- 
pendence of the amplitude Uo{t) can be determined by pro- 
jecting on the original velocity, i.e. 

Uo{t) = 2{U{x, t) sin fc, • x), (20) 

where the angular brackets denote volume averaging. The 
effective viscosity is connected with the decay rate A via 
A = i^tki, and A is measured as 

X=-{d\n\Uo{t)\/dt)t , (21) 
^ |http : //www ■ nordita ■ dk/ software /pencil- code| 
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Fig. 4. Time dependence of the vertical and radial velocity 
components at one point. Note that time is represented in vis- 
cous units, where 27r/fci is the vertical extent of the compu- 
tational domain. Evidently, the decay of the vertical epicyclic 
oscillations in is very small (lower panel), and there is no 
decay of Ux (upper panel). 

where {.)t denotes a time average over a suitable time interval 
when the decay is exponential. In forced turbulence it was 
found that vt ~ (0.8 . . . 0.9) x u^s/fef (Yousef et al. 2003). 

However, it turned out that in the present case the initial 
perturbations drive epicyclic oscillations that, once excited, 
do not easily decay. In Fig.l^we show the evolution of the ra- 
dial and vertical velocity components at one point (ux and Uz, 
respectively) in viscous units. Here the initial velocity per- 
turbation was chosen to be it ~ (0.1,0,0.1). The initial x 
velocity is compatible with the shearing periodic boundary 
conditions, but a uniform vertical velocity is not. This turns 
out to be not a problem, because viscosity is able to damp the 
initially produced sharp gradients that occur near the bound- 
ary to satisfy Uz = 0. 

Looking at Fig. E]it is clear that there is hardly any de- 
cay. On the other hand, if no perturbation is applied initially, 
some oscillations are still being generated by the turbulence, 
as can be seen from the temporal Fourier spectra of the ra- 
dial and vertical velocity components at one point in the do- 
main, |u2,(a;)p and \uz{uj)\'^, respectively. These are shown 
in Fig.lsjand Fig.|6lfor the cases with and without initial per- 
turbations that would initialize epicyclic oscillations. 

Indeed, in Fig.|5]one sees quite clearly the radial epicyclic 
frequency k — 0.4 in the radial velocity component and the 
vertical epicyclic frequency ( ~ 0.6 in the vertical velocity 
component. However, the latter seems to be shifted toward a 
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Fig. 5. Power spectra of Ux and Uz- The orbit of the initial 
state is perturbed, using u — (0.1,0,0.1), such that radial 
and vertical epicyclic motions are excited with frequencies 
K — n = 0.4 and ^ = 0.6 (marked by arrows in the two 
diagrams). The frequency of standing sound waves (p-modes) 
is labeled as ujp. The forcing amplitude is /q = 0.01. 
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Fig. 6. Same as Fig.|5j but without initial velocity perturba- 
tion corresponding to a perfectly circular orbit orbit. Note that 
the lower beat frequency, oji = — k = 0.2, appears to be 
excited in both horizontal and vertical velocity components, 
while the upper beat frequency, uj^ = C + ^ = appears 
only in the radial velocity component. 

somewhat higher frequency (about 0.7), which might be due 
to the interaction with the vertical p-mode frequency, ujp = 
CsTt/Lz « 0.78. 

It turns out that even when the oscillations are not present 
initially, some discrete frequencies are still being excited that 
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Fig. 7. Same as Fig.|5] but with magnetic field and no forcing, 
and in double-logarithmic representation. Thus, turbulence is 
generated by the MRJ. As in Fig.|5] an initial perturbation of 
amplitude 0. 1 has been applied. Note that for low frequencies 
the spectrum is nearly flat and lacks any discrete frequencies. 




lie near the lower and upper beat frequencies (uj\ — ( ~ k — 
0.2 and lu^ — C + k = 1; see Fig.|6l However, the match is 
only approximate, suggesting that the cause of the additional 
peaks in the spectra might not necessarily be related to beat 
phenomena. 

It is plausible that epicyclic oscillations in discs can be 
excited stochastically. This is analogous to the excitations of 
p-modes in the sun and in stars (Goldreich & Kumar 1990). 
Similar arguments may also be applied to discs in order to ex- 
plain the quasiperiodic oscillations (QPOs) in terms of reso- 
nances between different vertical and horizontal epicyclic fre- 
quencies (Abramowicz et al. 2003 a,b, Kato 2004, Kluzniaket 
al. 2004, Lee et al. 2004). 

Next we compare with a case where magnetic fields are 
included so that turbulence can be produced as a result of the 
MRI. No forcing function is therefore applied in the momen- 
tum equation. It turns out that in this case the power spectrum 
is more nearly flat and does not show discrete frequencies 
as in the purely hydrodynamic cases with an explicit forcing 
function; see Fig.0 

The absence of discrete frequencies in the velocity power 
spectrum of MRI-driven turbulence is surprising. In order to 
inspect further what happens after the initial kick that was 
imposed via the initial condition, we plot in Fig.|8]the time 
evolution for all three runs during the first 300 time units. 
(The total duration of these runs is around 10,000 time units.) 
It is clear from the figure that in the case of forced turbu- 
lence the amplitude of the epicycles remains dominant. In the 
case of MRI-driven turbulence, the epicycles are essentially 
swamped by the comparatively strong level of turbulence. By 
comparison, in the nonmagnetic case of forced turbulence 



Fig. 8. Comparison of the time series for the runs shown in 
Figs HJEl Only the beginning of the time series are shown. 
Note that in the last panel the MRI driven turbulence quickly 
wipes out the initial perturbation and no epicyclic oscillations 
are visible. 
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Fig. 9. Same as Fig. El but k = 0.22, ( = 0.33 « f k, and 
n = 0.4. 



some systematic oscillation pattern is visible even when there 
is no initial perturbation giving rise to epicycles. 
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Finally, we present a case where the shear parameter q is 
chosen to be 1.85, so the radial epicyclic frequency is now 
different from the rotation frequency, i.e. k = 0.22, while 
51 = 0.4. The vertical epicyclic frequency is chosen to be 
C = 0.33 ~ |k. It turns out that there is still a pronounced 
peak at 0.4 (— fl), while some of the earlier detected fre- 
quencies (0.2 and 1 .0), which are still present, lost their orig- 
inally anticipated interpretation; see Fig. |9] The vertical ve- 
locity has no longer a peak at ^ = 0.33, but instead at a 
higher frequency of about 0.5. However, it appears still pos- 
sible to interpret this as some interplay between C = 0.33 and 
ujp — 0.78. The frequency 0.2 is still visible in the spectra, 
but with the new values of k and ( it is no longer possible to 
explain this as the lower beat frequency between the two. Al- 
though this result is perhaps somewhat disappointing in terms 
of interpreting QPOs, these experiments highlight the general 
usefulness of using the shearing sheet approximation which 
allows these kinds of experiments to be carried out without 
difficulties. 

The above discussion is interesting in its own right and 
deserves certainly more attention. However, this method is 
obviously not well suited for determining turbulent transport 
coefficients. Therefore we now describe another method that 
has recently been developed in the hydromagnetic context for 
determining dynamo parameters. 

5. Tensorial turbulent resistivity 



The theory of turbulent resistivity is in many ways more de- 
veloped than the theory of turbulent viscosity. In dynamo the- 
ory it has recently become possible to determine quite accu- 
rately not only the turbulent resistivity, but also its full ten- 
sorial form and other components that can be non-dissipative 
and hence important for dynamo action. While in the hydro- 
dynamic case one is interested in the correlation UiUj, one is 
here interested in the correlation ujbj, or more specifically 
in the electromotive force £i — EijkUjbj. Assuming that the 
mean field is spatially smooth (which may not be the case in 
practice) one can truncate the expression for £i in terms of 
Bj and its derivatives after the first derivative, so one has 

Si = ctijBj + rjijkBj^k- (22) 

The components of tensor are usually quite easily de- 
termined from simulations by imposing a uniform magnetic 
field Bj and measuring the resulting electromotive force £i, 
so that aij = £i/Bj is obtained straightforwardly. The rea- 
son this works is because for a uniform field all derivatives 
of Bj vanish, so there are no higher order terms. Calcu- 
lating the components of rjijk is usually harder, especially 
when the mean field may no longer be smooth and its deriva- 
tives may vanish in places. A method that has been used for 
accretion disc turbulence is based on a fitting procedure of 
the measured mean field and the mean electromotive force 
to Eq. i22\ by calculating moments of the form {£iBj), 
{£iBj^k), as well as {BiBj) and (BiBj^k) (Brandenburg & 
Sokoloff 2002, Kowal et al. 2005). 

A general procedure for determining the full and rjijk 
tensors from a simulation is to calculate the electromotive 



force after applying test fields of different directions and 
with different gradients (Schrinner et al. 2005). In the fol- 
lowing we adopt xy averages, so the resulting mean fields 
depend only on z and t, and only B^ and By are non-trivial 
(Bz — because of the solenoidality B). Therefore, only the 
four components of and the four components of r/y 3 with 
i,j ~ 1,2 are non-trivial. Here, the numbers 1,2,3 refer to 
the cartesian x, y, z components. 

In the present case of one-dimensional mean fields it is 
advantageous to rewrite Eq. i22\ in the form 
£t = aijBj - fjijjj, i,j = 1,2, (23) 
where J = V x i? is the mean current density, and 
m = Vijk^jki (24) 
is the resistivity tensor operating only on the mean current 
density. In the special case of one-dimensional averages there 
is no extra information contained in the symmetric part of the 
Bj^k tensor that is not already contained in the components 
of J. In fact, the four components of rjij^ map uniquely to 
those of fjii via 

'/123 —'7113 
'72 2 3 —'7213 

This fact was also used in Brandenburg & Sokoloff (2002). 
The diagonal components of ijij correspond to turbulent re- 
sistivity, while its off-diagonal components can be responsi- 
ble for driving dynamo action [QxJ and W xJ effects; see, 
Radler (1969) and Rogachevskii & Kleeorin (2003, 2004), 
Radler & Stepanov (2005)]. Conversely, the diagonal com- 
ponents of the a tensor can be responsible for dynamo ac- 
tion while the off-diagonal components are responsible for 
non-regenerative turbulent pumping effects (Krause & Radler 
1980). It should be noted, however, that for linear shear flows 
Riidiger & Kitchatinov (2005) find that the signs of the rel- 
evant coefficients of fjij are such that dynamo action is not 
possible for small magnetic Prandtl numbers. 

In summary, in the present case of one-dimensional mean 
fields, B = B{z,t), there are altogether 4 + 4 unknowns. 
The idea is to calculate the electromotive force 

(26) 



( mi 


'712 ^ 






'722 / 





(25) 



t"^"^ =ux 



for the excess magnetic fluctuations, that are due to a 

given test field B^^''^\ where the labels p and q characterize 
the test field (p gives its nonvanishing component and q = 1 
or 2 stands for cosine or sine-like test fields. The calculation 
of the electromotive force requires solving simultaneously a 
set of equations of the form 

- \u + u)xb'^'-'^^ 



dt 



X 



2u{p,q) 



G (27) 



for each test field B^^ '^'' 



Here, G = V x [m x b^P''') - 
u X 

Ijip,'])] 

is a term nonlinear in the fluctuation. This term 
would be ignored in the first order smoothing approximation, 
but it can be kept in a simulation if desired. (In the present 
considerations it is neglected.) 

The four test fields considered in the present problem of 
one-dimensional mean fields are 

cos fciz \ / sin kiz 

1 , B*' " ={ I , (28) 



^(1.1) 
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B 



(2,1) 




B 



(2,2) 




As an example, consider the x component of £ 
and both values of q, 



(29) 



forp = 1 



-?(i,i) 
^1 



ail cos kiz — 77113 sin kiz, 



ti = ail sm fciz + 7^113 cos fciz. 



(30) 
(31) 



For p = 2, and/or for i = 2, one obtains a similar pair of 
equations with the same arrangement of cosine and sine func- 
tions. So, for each of the four combinations of i and j (= p) 



the set of two coefficient, and rjij^, is obtained as 




(32) 



where the matrix 

cos kiz 
sin kiz 



M 



(33) 



- sin kiz 
cos kiz 

is the same for each value of p and each of the two compo- 
nents i = 1, 2 of s'f''^'^ . Finally, fj is calculated using Eq. i24l . 
Note that det M = 1, so the inversion procedure is well be- 
haved and even trivial. 

The test field algorithm described above has been im- 
plemented in the PENCIL Code. The results are shown in 
Fig-Elfoi" the case of forced turbulence and in Fig.[^for the 
case of MRI-driven turbulence. We recall that in both cases 
the effect of stratification is included. 

The following main results can be summarized from this 
analysis. First, for forced turbulence the a effect (especially 
the yy component that is relevant for acting on the already 
strong toroidal By field) is positive in the northern hemi- 
sphere z > [as expected for cyclonic and anti-cyclonic 
events; see Parker (1979) and Krause & Radler (1980)]. How- 
ever, in the case of MRI-driven turbulence the sign of ayy is 
reversed, confirming the early results of Brandenburg et al. 
(1995). This may be explained in terms of strong flux tubes 
that are buoyant and therefore rising, but that also exhibit a 
converging flow from the B ■ 'VB tension force that points 
toward to strongest parts of the tube (Brandenburg 1998). The 
same result has been obtained by Riidiger & Pipin (2000) for 
magnetically driven turbulence using first order smoothing. 
Their result is specifically due to the dominance of the current 
helicity term and is not connected with the buoyancy term 
(Blackman 2005). 

Second, the off-diagonal components of the a tensor are 
such that they correspond to a turbulent pumping effect, 
£ = ...+7 X B, where 72 = ^(0^2; —a Theoretically the 
pumping velocity is expected to be in the direction of negative 
turbulent intensity (Roberts & Soward 1975). In the present 
case of forced turbulence the rms velocity is approximately 
independent of height, but the density decreases outwards, 
causing therefore an outward turbulent pumping effect (sec- 
ond panel of Fig.llO>. For MRI-driven turbulence, the density 
also deceases outward, but the rms velocity increases as the 
density decreases such that pit^ is approximately constant. 
This may be the reason for not having much of a vertical tur- 
bulent pumping effect in MRI-driven turbulence. 
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Fig. 10. Magnetic turbulent transport coefficients for forced 
turbulence. Solid and dashed lines refer respectively to the 
first and second quantity denoted on the corresponding axis. 
The error bars have been obtained by calculating the maxi- 
mum departure over the three possibilities obtained by con- 
sidering only 1/3 of the full time series. 



Third, the two diagonal components of fjij are in both 
cases positive (i.e. indeed diffusive, which is non-trivial!) and 
the two components, r]xx and fiyy, are almost equally big. [We 
recall that in Brandenburg & Sokoloff (2002) it was found 
that fjyy (responsible for diffusion of Bx) was very small, but 
this result was perhaps not accurate.] 

Fourth, the signs of the off-diagonal components of rjij 
are here such that they would not correspond to a dynamo ef- 
fect of the form £ = ... + S x J, where Sz = \{fjxy — Vyx)- 
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Fig. 11. Same as Fig.[To| but for MRJ-driven turbulence. 



Growing solutions require that the product of shear (here 
S = — ffi) and Sz is negative (Brandenburg & Subrama- 
nian 2005a,b). However, since < for forced turbu- 
lence (Fig. not , only decaying solutions are possible. This 
is in agreement with predictions by Riidiger & Kitchatinov 
(2005). For MRl-driven turbulence dz is within error bars ei- 
ther compatible with zero or positive in a few places. Thus, 
the W X J effect is perhaps possible here. However, for the 
full 2x2 tensor, fjij, it is important to consider its tenso- 
rial nature. It turns out that a necessary condition for dynamo 
action is 



fjyxk1{fjxyk1 + S) > (for dynamo action), 



(34) 



so the sign of fiyx is now crucial; see Appendix B for de- 
tails. Even then the simulations would not suggest dynamo 



action from the W x J effect. In any case, it is important to 
consider simulations with larger magnetic Reynolds. So far 
there is only the example of solar-like shear flow turbulence 
of Brandenburg (2005) where the x J effect is a leading 
candidate for explaining the generation of large scale fields 
in the absence of helicity. 

6. Conclusions 

The present investigations have demonstrated a number of 
new and previously unknown properties both for forced and 
MRl-driven turbulence in the shearing box approximation. In 
many respects these two cases are quite different. First of all, 
the fact that for MRl-driven turbulence most of the dissipa- 
tion happens near the upper and lower boundaries of the disc 
(i.e. in the disc corona) is peculiar to this case and is not in 
general true. Indeed, for non-magnetic discs most of the dis- 
sipation is expected near the midplane where the density is 
largest. Without entering the discussion about the reality of 
non-magnetic turbulence in accretion discs (e.g. Balbus et al. 
1996), we note that under some circumstances (e.g. in proto- 
stellar discs where the electric conductivity is poor) the MRl 
is not likely to operate, so less efficient mechanisms such 
as the inflow into the disc during its formation and vertical 
shear (Urpin & Brandenburg 1998) cannot be excluded as 
possible agents facilitating turbulence. Also the possibility of 
nonlinear instabilities (Richard & Zahn 1999, Chagelishvili 
et al. 2003, Afshordi et al. 2005) should be mentioned. In any 
case, comparing magnetic and non-magnetic cases is impor- 
tant in order to assess the potential validity of general turbu- 
lence concepts that have mainly been studied under forced 
non-magnetic conditions. 

There are a number of alternative ways of determining 
the turbulent viscosity in discs. Three methods have akeady 
been compared in the context of MRJ-driven turbulence by 
determining the stress either explicitly, from the heating rate, 
or from the radial mass accretion rate (Brandenburg et al. 
1996a). A rather obvious alternative is to consider the decay 
of an initial wave-like perturbation and to measure the decay 
of this signal. Although this method gives sensible results in 
the context of non-shearing and non-rotating flows (Yousef et 
al. 2003), in the present case it cannot be used for this pur- 
pose, because epicyclic oscillations are being initialized that 
hardly decay. In fact, at least in the hydrodynamic case with 
not too strong forcing there is clear evidence that epicyclic 
oscillations can actually be excited stochastically-very much 
like the p-modes in the sun. This may be important for under- 
standing the quasi-periodic oscillations discovered recently 
in some pulsars (e.g., Lee et al. 2004), and in particular the 
connection with the possibility of exciting resonances at the 
beat frequencies for non-equal vertical and radial epicyclic 
frequencies. 

Given that the turbulence is in general non-isotropic (ow- 
ing to shear and rotation, as well as stratification), the turbu- 
lent transport coefficients are in general also non-isotropic. 
For practical calculations of toroidally averaged accretion 
flows (e.g., Kley et al. 1993, Igumenshchev et al. 1996) it 
is therefore essential to know the full tensorial structure as 
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well as other possibly non-diffusive contributions. Signifi- 
cant progress has been made in determining the full tensorial 
structure of the turbulent resistivity and alpha effect. Some- 
what surprisingly, it turned out that the two diagonal compo- 
nents of the resistivity tensor are nearly equal. Furthermore, 
in the non-MRI case the off-diagonal components can even 
given rise to dynamo action (the so-called W x J or shear- 
current effect). As far as the a effect is concerned, an earlier 
result about the different signs for MRI and forced turbulence 
is confirmed. 

Similar investigations can probably also be carried out 
for determining the tensorial nature of turbulent viscosity 
and possibly other non-diffusive effects that are known in 
other circumstances (e.g., Riidiger & Hollerbach 2004). Most 
important is perhaps the investigation of stochastically ex- 
cited epicyclic oscillations in connection with the kilohertz 
quasiperiodic oscillations found in some pulsars. Obviously, 
more systematic investigations should be carried out to deter- 
mine the dependence on forcing amplitude of the turbulence 
and to check whether similar oscillations are also possible for 
MRI-driven turbulence. 
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Appendix A: The forcing function 

Ror completeness we specify here the forcing function used in the 
present paper^ . It is defined as 

f{x, t) = Re{Ar/j^(^j exp[ife(t) ■ x + i0(t)]}, (Al) 

where x is the position vector. The wavevector k{t) and the ran- 
dom phase — TT < (f){t) < n change at every time step, so f{x,t) 
is (5-correlated in time. Ror the time-integrated forcing function to 
be independent of the length of the time step St, the normalization 
factor A'^ has to be proportional to St^^^'^ . On dimensional grounds 

' This forcing function was also used by Brandenburg (2001), but 
in his Eq. (5) the factor 2 in the denominator should have been re- 
placed by \/2 for a proper normalization. 
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it is chosen to be A'^ = foCs{\k\cs/5ty^^ , wliere fo is a nondi- 
mensional forcing amplitude. Tlie value of the coefficient fo is cho- 
sen such that the maximum Mach number stays below about 0.5; in 
practice this means fo = 0.01 . . . 0.05, depending on the average 
forcing wavenumber. At each timestep we select randomly one of 
many possible wavevectors in a certain range around a given forc- 
ing wavenumber. The average wavenumber is referred to as kf. Two 
different wavenumber intervals are considered: 1...2 for kf = 1.5 
and 4. 5. ..5. 5 for kf — 5. We force the system with transverse heli- 
cal waves, 

fk = {kx e) /Vfe2-(fc-e)2, (A2) 
where e is an arbitrary unit vector not aligned with fe; note that 

Appendix B: Dispersion relation for tensorial fj 

The dispersion relation for the dynamo problem with 5 effect is 
(Brandenburg & Subramanian 2005) 

A± = -r/Tfc^ ± V-(fc- 5)2fc2 - {k-S)Sk„ (Bl) 

where A is the growth rate, k^ and fez are the wavenumbers in the 
X and z directions, fc^ = ki + k1, and rjT = rj + rjt the sum of 
microscopic and turbulent diffusion. 

We now consider the one-dimensional problem (fea: — 0) and 
treat all four components of fjij as independent, so we have the fol- 
lowing set of two partial differential equations, 

= {ri + fjyy)B'^ ~ fjyxBy, (B2) 

'iy ^ {ri + flxx)By ~fj:,yB" + SBx, (B3) 

where S = — |f2 is the shear coefficient. Using the abbreviation 
""-'j ~ iv^ij + Vij)^^' this leads to the dispersion relation 

A± = —^{nxx + riyy) ± j(nxx — nyy)'^ + nyx(nxy + S), (B4) 

which reduces to the one-dimensional form of Eq. <B1> when fixx ~ 
Vyy = Vt and fixy = —Vyx = S. 
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